library(ggplot2)
library(patchwork)
'%!in%' <- function(x,y)!('%in%'(x,y))

# set to desired working directory
# setwd("")

# Start log file for script -----------------------------------------------
TeachingDemos::txtStart("./log files/figureA22.txt")

# Load Data ---------------------------------------------------------
# individual data
load("./data/impdat20.Race.WeekWeight.rdata")
load("./data/impdat20.Age.WeekWeight.rdata")
load("./data/impdat20.Arab.WeekWeight.rdata")
load("./data/impdat20.Asian.WeekWeight.rdata")
load("./data/impdat20.Disab.WeekWeight.rdata")
load("./data/impdat20.GendCar.WeekWeight.rdata")
load("./data/impdat20.GendSci.WeekWeight.rdata")
load("./data/impdat20.NatAm.WeekWeight.rdata")
load("./data/impdat20.Pres.WeekWeight.rdata")
load("./data/impdat20.Religion.WeekWeight.rdata")
load("./data/impdat20.Sexuality.WeekWeight.rdata")
load("./data/impdat20.Weapons.WeekWeight.rdata")
load("./data/impdat20.Weight.WeekWeight.rdata")

# Race
impdat20.WeekWeight$race_black <- ifelse(impdat20.WeekWeight$race5 == "Black", 1, 0)
impdat20.WeekWeight$race_black[which(is.na(impdat20.WeekWeight$race5))] <- NA

impdat20.WeekWeight$race_black <- ifelse(impdat20.WeekWeight$race5 == "Black", 1, 0)
impdat20.WeekWeight$race_black[which(is.na(impdat20.WeekWeight$race5))] <- NA

impdat20.Age.WeekWeight$race_black <- ifelse(impdat20.Age.WeekWeight$race5 == "Black", 1, 0)
impdat20.Arab.WeekWeight$race_black <- ifelse(impdat20.Arab.WeekWeight$race5 == "Black", 1, 0)
impdat20.Asian.WeekWeight$race_black <- ifelse(impdat20.Asian.WeekWeight$race5 == "Black", 1, 0)
impdat20.Disab.WeekWeight$race_black <- ifelse(impdat20.Disab.WeekWeight$race5 == "Black", 1, 0)
impdat20.GendCar.WeekWeight$race_black <- ifelse(impdat20.GendCar.WeekWeight$race5 == "Black", 1, 0)
impdat20.GendSci.WeekWeight$race_black <- ifelse(impdat20.GendSci.WeekWeight$race5 == "Black", 1, 0)
impdat20.NatAm.WeekWeight$race_black <- ifelse(impdat20.NatAm.WeekWeight$race5 == "Black", 1, 0)
impdat20.Pres.WeekWeight$race_black <- ifelse(impdat20.Pres.WeekWeight$race5 == "Black", 1, 0)
impdat20.Religion.WeekWeight$race_black <- ifelse(impdat20.Religion.WeekWeight$race5 == "Black", 1, 0)
impdat20.Sexuality.WeekWeight$race_black <- ifelse(impdat20.Sexuality.WeekWeight$race5 == "Black", 1, 0)
impdat20.Weapons.WeekWeight$race_black <- ifelse(impdat20.Weapons.WeekWeight$race5 == "Black", 1, 0)
impdat20.Weight.WeekWeight$race_black <- ifelse(impdat20.Weight.WeekWeight$race5 == "Black", 1, 0)
impdat20.Age.WeekWeight$race_black[which(is.na(impdat20.Age.WeekWeight$race5))] <- NA
impdat20.Arab.WeekWeight$race_black[which(is.na(impdat20.Arab.WeekWeight$race5))] <- NA
impdat20.Asian.WeekWeight$race_black[which(is.na(impdat20.Asian.WeekWeight$race5))] <- NA
impdat20.Disab.WeekWeight$race_black[which(is.na(impdat20.Disab.WeekWeight$race5))] <- NA
impdat20.GendCar.WeekWeight$race_black[which(is.na(impdat20.GendCar.WeekWeight$race5))] <- NA
impdat20.GendSci.WeekWeight$race_black[which(is.na(impdat20.GendSci.WeekWeight$race5))] <- NA
impdat20.NatAm.WeekWeight$race_black[which(is.na(impdat20.NatAm.WeekWeight$race5))] <- NA
impdat20.Pres.WeekWeight$race_black[which(is.na(impdat20.Pres.WeekWeight$race5))] <- NA
impdat20.Religion.WeekWeight$race_black[which(is.na(impdat20.Religion.WeekWeight$race5))] <- NA
impdat20.Sexuality.WeekWeight$race_black[which(is.na(impdat20.Sexuality.WeekWeight$race5))] <- NA
impdat20.Weapons.WeekWeight$race_black[which(is.na(impdat20.Weapons.WeekWeight$race5))] <- NA
impdat20.Weight.WeekWeight$race_black[which(is.na(impdat20.Weight.WeekWeight$race5))] <- NA

# * Age ---------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Young_Good_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Age.WeekWeight, 
                    subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25",
                    weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Age"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
age_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# * Arab ------------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.OtherPeople_Good_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Arab.WeekWeight, 
                    subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", 
                    weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Arab"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
arab_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# * Asian -----------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.White_American_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Asian.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Asian"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
asian_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Disability --------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Abled_Good_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Disab.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Disability"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
disability_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# * Gender-Career ---------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Male_Career_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.GendCar.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Gender-Career"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
gendcar_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# Gender-Science ----------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Male_Science_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.GendSci.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Gender-Science"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
gendsci_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Native American ---------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.White_American_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.NatAm.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Native American"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
natam_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# President ---------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Trump_Good_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Pres.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "President"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
pres_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Religion ----------------------------------------------------------------

# ^ Christianity-Islam ----------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Christianity_Ci_Good_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Religion.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Christianity-Islam"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
relig_ci_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# ^ Christianity-Judaism --------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Christianity_CJ_Good_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Religion.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Christianity-Judaism"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
relig_cj_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())



# ^ Judaism-Islam ---------------------------------------------------------

# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Judaism_Ji_Good_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Religion.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Judaism-Islam"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
relig_ji_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Sexuality ---------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Straight_Good_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Sexuality.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Sexuality"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
sexuality_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Weapons -----------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Black_Weapons_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Weapons.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Weapons"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
weapons_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Weight -----------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_bw_dummy <- lm(D_biep.Thin_Good_all ~ 
                      as.factor(week)*race_black + 
                      ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Weight.WeekWeight, subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_bw_dummy)

# ** Figure ----------------------------------------------------------
title <- "Weight"
weeks <- names(table(m1_bw_dummy$model$`as.factor(week)`))

# Race
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_wht <- predict(m1_bw_dummy, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw_dummy, pdat_blk, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01"

# Plotting
weight_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# Race IAT for Placebo Comparison -----------------------------------------
m1_bw <- lm(D_biep.White_Good_all ~ 
              as.factor(week)*race_black + ideo_mod + ideo_con +
              age_cat4 + edu_cat4 + sex + 
              census_region + broughtwebsite2,
            data = impdat20.WeekWeight, weights = WeekWeights,
            subset = (race5 == "White" | race5 == "Black") & date_mdy != "2020-05-25")
summary(m1_bw)

weeks <- names(table(m1_bw$model$`as.factor(week)`))

# Data Frames for Predictions
pdat_wht <- data.frame(week = weeks,
                       race_black = 0,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_blk <- data.frame(week = weeks,
                       race_black = 1,
                       ideo_mod = 0,
                       ideo_con = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

weeks <- names(table(m1_bw$model$`as.factor(week)`))
preds_wht <- predict(m1_bw, pdat_wht, se.fit = T)
preds_blk <- predict(m1_bw, pdat_blk, se.fit = T)

# combine into DF
p_dat <- data.frame(week = weeks,
                    race = c(rep("White", length(weeks)),
                             rep("Black", length(weeks))),
                    pred = c(preds_wht$fit, preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit))
p_dat_fill <- data.frame(week = unique(impdat20.WeekWeight$week)[which(unique(impdat20.WeekWeight$week) %!in% weeks)],
                         race = c(rep("White", length(unique(impdat20.WeekWeight$week)[which(unique(impdat20.WeekWeight$week) %!in% weeks)])),
                                  rep("Black", length(unique(impdat20.WeekWeight$week)[which(unique(impdat20.WeekWeight$week) %!in% weeks)]))),
                         pred = NA,
                         se = NA)
p_dat <- rbind(p_dat, p_dat_fill)
p_dat$race <- factor(p_dat$race,
                     levels = c("Black",
                                "White"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1)] <- "2020-01-01" # from 00 weeks


# Plot
race_fig <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2, color = "grey") +
  geom_linerange(aes(ymin = pred - 1.65*se, 
                     ymax = pred + 1.65*se),
                 lwd = .75, color = "grey") +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred), se = FALSE, color = "black") +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred), se = FALSE, color = "black") +
  facet_grid(race~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "") +
  scale_y_continuous(labels = numform::ff_num(zero = 0, digits = 2)) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 20),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 20),
        strip.background = element_blank())

lhs_fit <- ggplot_build(race_fig)$data[[3]]
lhs_fit$y[which(lhs_fit$PANEL == 1)]
lhs_fit$y[which(lhs_fit$PANEL == 2)]

rhs_fit <- ggplot_build(race_fig)$data[[4]]
rhs_fit$y[which(rhs_fit$PANEL == 1)]
rhs_fit$y[which(rhs_fit$PANEL == 2)]


race_preds1 <- data.frame(iat = "race",
                          out = "loess",
                          Black = rhs_fit$y[which(rhs_fit$PANEL == 1)][1] - lhs_fit$y[which(lhs_fit$PANEL == 1)][length(lhs_fit$y[which(lhs_fit$PANEL == 1)])],
                          White = rhs_fit$y[which(rhs_fit$PANEL == 2)][1] - lhs_fit$y[which(lhs_fit$PANEL == 2)][length(lhs_fit$y[which(lhs_fit$PANEL == 2)])])
race_preds2 <- data.frame(iat = "race",
                          out = "predicted",
                          Black = p_dat[which(p_dat$week == "2020-05-25" & p_dat$race == "Black"), "pred"] - p_dat[which(p_dat$week == "2020-05-18" & p_dat$race == "Black"), "pred"],
                          White = p_dat[which(p_dat$week == "2020-05-25" & p_dat$race == "White"), "pred"] - p_dat[which(p_dat$week == "2020-05-18" & p_dat$race == "White"), "pred"])
race_preds <- rbind(race_preds1, race_preds2)
race_preds <- reshape::melt(race_preds, id = c("iat", "out"), variable_name = "race")


# Plot -----------------------------------------------------
weeks <- unique(p_dat$week)

obs_diffs <- as.data.frame(array(NA, c(14, 3)))
names(obs_diffs) <- c("iat", "Black", "White")

placeb_diffs <- as.data.frame(array(NA, c(1, 4)))
names(placeb_diffs) <- c("iat", "week", "Black", "White")

x_diffs <- as.data.frame(array(NA, c(length(weeks), 4)))
names(x_diffs) <- c("iat", "week", "Black", "White")

plot_names <- grep("_plot", ls(), value = T)
plot_names <- plot_names[which(plot_names != "facet_plot")]

# loop across plot objects to pull predicted D-scores for observed differences and pseudo differences
for(i in 1:length(plot_names)){
  # metadata storing
  x_diffs$iat[1] <- gsub("_plot", "", plot_names[i])
  x_diffs$week[1] <- weeks[1]
  
  # pull observed George Floyd Difference
  obs_diffs$iat[i] <- gsub("_plot", "", plot_names[i])
  
  lhs_fit <- ggplot_build(get(plot_names[i]))$data[[1]]
  rhs_fit <- ggplot_build(get(plot_names[i]))$data[[2]]
  
  obs_diffs$Black[i] <- rhs_fit$y[which(rhs_fit$PANEL == 1)][which(p_dat$week == "2020-05-25")[1]] - lhs_fit$y[which(lhs_fit$PANEL == 1)][which(p_dat$week == "2020-05-18")[1]]
  obs_diffs$White[i] <- rhs_fit$y[which(rhs_fit$PANEL == 2)][which(p_dat$week == "2020-05-25")[1]] - lhs_fit$y[which(lhs_fit$PANEL == 2)][which(p_dat$week == "2020-05-18")[1]]

  # other differences but for the treatment period
  for(j in 2:length(weeks)){
    x_diffs$iat[j] <- gsub("_plot", "", plot_names[i])
    x_diffs$week[j] <- weeks[j]
    
    if(weeks[j] <= "2020-05-18"){
      x_diffs$Black[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 1)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 1)][which(p_dat$week == weeks[j-1])[1]]
      x_diffs$White[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 2)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 2)][which(p_dat$week == weeks[j-1])[1]]
    }
    if(weeks[j] == "2020-05-25" | weeks[j] == "2020-06-01"){
      # empty placeholder
    }
    if(weeks[j] > "2020-06-01"){
      x_diffs$Black[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 1)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 1)][which(p_dat$week == weeks[j-1])[1]]
      x_diffs$White[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 2)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 2)][which(p_dat$week == weeks[j-1])[1]]
    }
  }
  
  # store results
  placeb_diffs <- rbind(placeb_diffs, x_diffs)
}

# Reformatting and cleaning up
placeb_diffs <- placeb_diffs[-1,]
placeb_diffs <- reshape::melt(placeb_diffs, id = c("iat", "week"), variable_name = "race")
obs_diffs <- reshape::melt(obs_diffs, id = c("iat"), variable_name = "race")

# Getting p-values 

# Data frame to store p-values
p_values <- obs_diffs

for(i in 1:nrow(p_values)){
  # store test statistic
  ob_val <- obs_diffs$value[i]
  # extract appropriate simulations for reference distribution
  ref_dist <- subset(placeb_diffs, iat == p_values$iat[i] & race == p_values$race[i],
                     select = value)
  ref_dist <- ref_dist$value[which(!is.na(ref_dist$value))]
  # store number of simulations
  n_sim <- length(ref_dist)
  # two-tailed p-value
  p_values$value[i] <- sum(abs(ref_dist) > abs(ob_val))/n_sim
}
p_values$x <- -.15
p_values$y <- 15
p_values$p <- paste0("p=", round(p_values$value, 2))

unique(placeb_diffs$iat)

labs <- c("Age", "Arab", "Asian", "Disability",
          "Gender-\nCareer", "Gender-\nScience",
          "Native\nAmerican", "President", 
          "Religion\nChristian-\nIslam", "Religion\nChristian-\nJudaism",
          "Religion\nJudaism-\nIslam", "Sexuality", "Weapons", "Weight")

placeb_diffs$iat <- factor(placeb_diffs$iat, labels = labs)
p_values$iat <- factor(p_values$iat, labels = labs)
obs_diffs$iat <- factor(obs_diffs$iat, labels = labs)


# Plotting
ggplot(placeb_diffs, aes(x = value)) +
  geom_histogram(fill = "grey") +
  geom_vline(data = obs_diffs, aes(xintercept = value)) +
  geom_text(data = p_values, aes(label = p, x = x, y = y)) +
  facet_grid(iat~race, scales = "free_x") +
  scale_x_continuous(labels = numform::ff_num(zero = 0, digits = 2), limits = c(-.45, .45)) +
  theme_bw() +
  labs(x = "Predicted IAT Score Difference", y = "Count") +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 16),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 16),
        strip.background = element_blank())
ggsave("./figures/FigureA22.pdf", height = 16, width = 8)

# End log file for script -------------------------------------------------
TeachingDemos::txtStop()